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Abstract 


The use of decohesion elements for the simulation of delamination in composite materials 
is reviewed. The test methods available to measure the interfacial fracture toughness 
used in the formulation of decohesion elements are described initially. After a brief 
presentation of the virtual crack closure technique, the technique most widely used to 
simulate delamination growth, the formulation of interfacial decohesion elements is 
described. Problems related with decohesion element constitutive equations, mixed-mode 
crack growth, element numerical integration and solution procedures are discussed. 
Based on these investigations, it is concluded that the use of interfacial decohesion 
elements is a promising technique that avoids the need for a pre-existing crack and pre- 
defined crack paths, and that these elements can be used to simulate both delamination 
onset and growth. 


Introduction 

The fracture process of high performance composite laminates is quite complex, involving both 
intralaminar damage mechanisms (e.g. matrix cracking, fiber fracture) and interlaminar damage 
(delamination). An example of a failure with interactive modes is illustrated in Figure 1 . Although some 
progress has been made lately in the development of accurate analytical tools for the prediction of 
intrala mi nar damage growth, similar tools for delamination are still not available, and thus delamination is 
generally not considered in damage growth analyses. Without the delamination failure mode, the 
predictive capabilities of progressive failure analyses will remain limited. 



Figure 1. Interaction between intralaminar and interlaminar damage mechanisms [1], 
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Figure 2. Experiment illustrating stiffener-flange debonding. 


Delamination is one of the predominant forms of failure in laminated composites due to the lack of 
reinforcement in the thickness direction. Delamination as a result of impact or a manufacturing defect can 
cause a significant reduction in the compressive load-carrying capacity of a structure. The stress gradients 
that occur near geometric discontinuities such as ply drop-offs, stiffener terminations, s ki n-stiffener 
flange interfaces (Figure 2), bonded and bolted joints, and access holes promote delamination initiation, 
trigger intraply damage mechanisms, and can result in a significant loss of structural integrity. 

The analysis of delamination is commonly divided into the study of the initiation and the analysis of 
the propagation of an already initiated area. Delamination initiation analysis is usually based on stresses 
and use of criteria such as the quadratic interaction of the interlaminar stresses in conjunction with a 
characteristic distance [2,3]. This distance is a function of specimen geometry and material properties, 
and its determination always requires extensive testing. 

Crack propagation, on the other hand, is usually predicted using the Fracture Mechanics approach. The 
Fracture Mechanics approach avoids the difficulties associated with the stress singularity at a crack front 
but requires the presence of a pre-existing delamination whose exact location may be difficult to 
determine. When used in isolation, neither the strength-based approach nor the Fracture Mechanics 
approach is adequate for a comprehensive analysis of progressive delamination failure. 

The objective of this report is to examine the several aspects of simulating the delamination, including 
a method that combines elements of the strength and Fracture Mechanics approaches. A description of 
test methods that have been proposed to obtain the required material properties is presented. The use of 
the virtual crack closure technique and interfacial decohesion elements for the simulation of delamination 
failure is discussed. The main problems associated with the formulation and implementation of interfacial 
decohesion elements in a finite element model are presented. 

Measurement of Interfacial Fracture Toughness 

It is essential for developers of computational methods to initiate the analysis with a careful 
observation of the process to be simulated. In the particular case of delamination growth, it is considered 
that experimental observation of the physical mechanisms occurring at the interfaces during crack growth 
provides essential information for the simulation of the process. 

The different test methods used for the determination of interlaminar fracture toughness provide an 
effective means to observe the damaged surfaces under different loading conditions. A careful 
examination of the characteristics of each test method is required to define the appropriate information to 
be used in the numerical models. 
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Mode I 


The geometry used to determine the interlaminar fracture toughness in Mode I (G IC ) is the double- 
cantilever beam (DCB) specimen. This specimen is made of a unidirectional fiber-reinforced laminate 
containing a thin insert at the mid-plane near the loaded end (Figure 3). 


Applied force 


V- 

Insert 


Figure 3. DCB test specimen. 


Davies et al. [4] suggested that the major problems associated with this test are the occurrence of fiber 
bridging across the crack as the crack moves above and below bundles of fibers and the fact that the R- 
curves associated with this phenomenon are not intrinsic material properties, but they frequently depend 
on specimen stiffness. However, OBrien [5] noted that fiber bridging, and consequently the R-curve 
effect shown in Figure 4, are artifacts of the unidirectional DCB specimen and do not occur in structural 
composite laminates. 

The most appropriate load introduction mechanism, the film thickness of the insert and the data 
reduction methods for calculating G/c values were proposed as a result of a comprehensive round-robin 
test program [6]. Although the standard test method limits the DCB specimens to unidirectional (0 s ) 
laminates, it is expected that delaminations in composite structures preferentially arise on interfaces 
between layers with different orientations. Hence, it is important to have experimental information 
concerning interlaminar fracture toughness for interfaces between other than 0 s plies. Preliminary results 
have shown that unidirectional composites provide the most conservative values for Gic [4], 

The main problem associated with performing DCB tests on interfaces other than 0 S /0 S is that the 
delamination branches to interfaces away from the midplane [6-8], Experimental results [7, 8] have 
shown higher values of G/c for interfaces other than 0 S /0 S . However, in all of the tests in Refs. 7 and 8 
crack jumping was observed and the data was analyzed using standard procedures. This phenomenon has 
three consequences: firstly, the fracture energy reported is affected by the intralaminar fracture toughness 
due to the intralaminar matrix cracks. Secondly, as the crack moves from one interface to another, the 
stiffness of the laminates above and below the crack changes, and this effect should be taken into account 
in the expressions used to determine G IC . Finally, a DCB test using arms with unequal bending stiffness 
will also promote Mode II loading [9], 

The modified DCB, a new test method designed to overcome these problems, has been proposed by 
Robinson and Song [10], Investigations using [-45 s /0 s /(45 s )2/0 e /-45 s /45 s /0 s /(-45 s )2/0 2 /45 s )] s , 
[-45 2 /0 2 /(45 2 )2/0 2 /-45 2 /45 s /0 a /(-45 2 )2/0 2 /45 a )] 2 and [ -45 s /45 s /(0 e )2/45 2 /-45 2 /4 5 2 /-4 5 2 /(0 e )2/-45 2 /45 2 )] 2 edge 
delaminated specimens to determine fracture toughness at +45 2 /-45 2 and +45 2 /+45 2 interfaces have shown 
a suppression of the crack jumping that occurs in the conventional DCB specimen. Preliminary tests [10] 
have also resulted in significantly higher values of G [C for multidirectional laminates compared with 
unidirectional laminates. 
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Figure 4. R-curve effect on a DCB test of a (0 e )24 CFRP laminate [7]. 


Mode II 

There are four test specimens being considered for the measurement of interlaminar fracture 
toughnesses under Mode II loading [4, 5]: the end notched flexure specimen (ENF), the stabilized end 
notched flexure specimen (SENF), the four point end notched flexure specimen (4ENF), and the end 
loaded split specimen (ELS). Although the ENF test specimen (Figure 5) has received the most attention, 
it has problems associated with the unstable crack propagation for short crack lengths. One way to 
stabilize the ENF test is by feedback load control of the test machine [4, 11]. 


Applied load 




z: 


t 


Insert 


t 


Figure 5. ENF test specimen. 


According to Davies [4], the main problems associated with Mode II testing are the definition of the 
type of starter defect, the definition of crack initiation, the stability of the test, frictional effects on the 
crack faces and the data analysis. Furthermore, the Mode II fracture toughness is typically much higher 
and has more scatter than the Mode I fracture toughness [12]: for epoxy matrix composites the GucIGic 
ratios typically exceed a value of 2. 

Observing the fracture surfaces that formed under Mode I and Mode II, O'Brien [12] concluded that 
under Mode II loading the fracture surface exhibits a rough fracture plane, including ‘hackles’ (Figure 6), 
whereas the Mode I fracture surface has a clean cleavage plane. O'Brien also concluded that for both 
composite laminates and structural elements, failure occurs at a location where Mode I accounts for at 
least half of the total G at failure, so Mode I and mixed-mode interlaminar fracture toughnesses are the 
most relevant for predicting delamination failure. 
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Figure 6. Microcracks forming the shear fracture surface [13]. 


Mode III 

The majority of the research on the determination of interlaminar fracture toughness has been 
performed for Mode I and Mode II loading. For a complete characterization of the fracture process, Mode 
III delamination tests are also required. 


I Applied load 


I -^-1 

f Insert t 

Figure 7. MMB test specimen. 


A test that has received much interest is the edge crack torsion (ECT), which is based on the out-of- 
plane torsion of a cracked (90 2 /(±45 2 )n/( + 45 2 )n/90 2 ) s specimen [14]. According to Davies [4], the same 
concerns occurring in Mode II loading regarding type of defect, initiation definition, and friction need to 
be addressed for the ECT specimen. Furthermore, the twist-bending coupling due to the presence of off- 
axis plies leads to a Mode II contribution to the total strain energy release rate that needs to be identified 
[15], Using a shear deformation theory, Li and Wang [16] have shown that the Mode II contribution to the 
total strain energy release rate was negligible. 

The influence of the transverse shear modulus G 23 on the Mode III toughness was investigated by Li 
and O'Brien [15], G 23 is required for the analysis [16] and it is difficult to determine experimentally. It 
was concluded that a higher value of G23 yields lower Mode III toughness estimations and that assuming 
G 2 jto be equal to Gn results in a conservative estimate for the initial delamination length. 

Mixed-Mode I and II 

The most widely used specimen for mixed-mode fracture is the mixed-mode bending (MMB) specimen, 
shown in Figure 7, which was proposed by Reeder and Crews [17], This specimen was later re-designed 
to minimize geometric nonlinearities [18], 

According to Davies [4], the main advantages of the MMB test method are the possibility of using 
virtually the same specimen geometry for Mode I tests and varying the mixed mode ratio from pure Mode 
I to pure Mode II. 
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Numerical Simulation of Delamination 


Virtual Crack Closure Technique 

The virtual crack closure technique (VCCT) proposed by Rybicki and Kanninen [19] has been used to 
predict delamination growth in composite materials. This technique is based on Irwin’s assumption that 
when a crack extends by a small amount, the energy absorbed in the process is equal to the work required 
to close the crack to its original length. The energy release rates can then be computed from the nodal 
forces and displacements obtained from a finite element model. 

The approach can be computationally effective when sufficiently refined meshes are used and when all 
the elements at the crack tip have the same dimensions in the crack growth direction. Under these 
conditions, the energy release rates can be obtained from only one analysis. 

Koning et al. [20] have used the VCCT to predict delamination growth in plates containing a circular 
delamination of 10-mm diameter and loaded in tension and compression. A three-dimensional layered 
element with eight nodes was used and the Mode I, Mode II, and Mode III energy release rates were 
computed along the delamination front. The location of maximum predicted energy release rates, and 
hence the delamination growth, were in good agreement with experimental results. 

The use of three-dimensional elements to determine energy release rates using the VCCT can 
overcome the dependence of the results on element shear deformation assumptions. However, these 
elements lead to computationally intensive models. In order to minimize this problem, Krueger [21] has 
proposed combining shell and three-dimensional elements. In this approach, the solid three-dimensional 
elements are used in the immediate vicinity of the delamination front. The connection between the shell 
and the three-dimensional elements was performed by enforcing the appropriate translations and rotations 
at the interface. The procedure was used in the simulation of DCB, ENF and single leg bending (SLB) 
specimens for delaminations located at 0 2 /0 2 and +3() 2 /-30 2 interfaces. 

The results of the DCB specimen test clearly show a higher energy release rate at the center of the 
0 2 /0 2 specimen compared to the specimen edges due to the anticlastic bending effect. This effect was even 
more pronounced in the specimens containing +30 2 /-30 2 interfaces due to the lower bending stiffness of 
its arms. 

The Gn values for the ENF specimens are constant across most of the specimen width with higher 
values in the vicinity of the edges, accompanied by a localized Mode III contribution. The accuracy of the 
analysis depends on the size of the region modeled with three-dimensional elements: results similar to the 
ones from a full three-dimensional analysis were obtained with the local zone extended to three times the 
specimen thickness [21], Furthermore, for refined meshes, 8-node brick elements combined with 4-node 
shell elements provided identical results to a model with 20-node brick elements combined with 8-node 
shell elements. However, care must be taken when simulating the specimen with +30 2 /-30 2 interfaces, 
since experimental results obtained for interfaces other than 0 2 /0 2 have shown delamination jumping from 
one interface to another [7], and the procedure proposed is unable to take this effect into consideration. 

The VCCT has also been used to determine energy release rates for delaminations originating from 
matrix cracks in skin-stringer interfaces when subjected to arbitrary load conditions [22], Using 8-node 
plane strain elements with reduced integration, a G versus delamination length a relation was obtained. 
From these results, it was possible to obtain information concerning the stability of delamination growth. 

Care must be taken in selecting the element size at the crack tip when using the VCCT to simulate 
delamination. Raju et al. [23] have shown that the individual components of the energy release rate does 
not converge when the ratio of the size of delamination tip element to the ply thickness decreases. This 
effect does not occur for the total energy release rate, and it is due to the oscillatory part of the stress 
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singularity occurring in cracks between dissimilar media. It was shown that when using delamination tip 
elements with a size of one-quarter to one-half of the ply thickness, the individual components obtained 
agreed well with the results from a model that had a resin-rich layer incorporated between the plies [23], 

Although providing valuable information concerning onset and stability of delamination growth, the 
use of VCCT to simulate dela mi nation growth requires complex moving mesh techniques [24] to advance 
the crack front when the local energy release rate exceeds a critical value. Furthermore, an initial 
delamination must be defined and, for certain geometries and load cases, the location of the delamination 
front might be difficult to determine. 

Interfacial Decohesion Elements 

One approach to overcome the above limitations of the VCCT is the use of interfacial decohesion 
elements placed between the composite material layers. Interfacial decohesion elements combine a stress 
based formulation with a Fracture Mechanics based formulation and are used to define the non-linear 
constitutive law of the material at the interface between laminae. 

The concept of decohesion elements is based on a Dugdale-Barenblatt type cohesive zone [25, 26] and 
was first applied to the analysis of concrete cracking by Hillerborg et al. [27]. The concept of cohesive 
zone models has also been used by Needleman [28] to describe the process of void nucleation from 
inclusions, to simulate fast crack growth in brittle solids [29], and for intersonic crack growth under shear 
loading [30], 

Needleman [28] considered cohesive zone models particularly attractive when interfacial strengths are 
relatively weak compared with the adjoining material, as in the case of composite laminates. The 
proposed constitutive equations for the interface are phenomenological mechanical relations between the 
tractions and interfacial separations such that, with increasing interfacial separation, the tractions across 
the interface reaches a maximum, decreases and vanishes when complete decohesion occurs, as illustrated 
in Figure 8. 






7 


\ / 






Figure 8. Normal and shear tractions across the cohesive surface as a function of relative 
displacements [29], 
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The decohesion response was specified in terms of a potential relating the interfacial tractions and the 
relative tangential and normal displacements across the interface [29], The resulting work of normal 
separation and tangential separation can be related to the critical values of energy release rates. 

In fact, the cohesive zone formulation is identical to Griffith’s theory of fracture: considering the 
integral J proposed by Rice [31]: 


J= \ (wdy-T.—ds) (1) 

Jr dx 

where r is any contour surrounding the notch tip, u is the displacement vector, T is the traction vector 
defined accordingly to the outward normal along /] and w is the strain energy density. Evaluating the 
integral J on the contour /"defined by the cohesive zone (Figure 9): 


1 O 

J = f ( wdy - a (S) — dx) 

Jr rlv 


dx 


( 2 ) 


Taking into account that dy=0 along the selected contour [31]: 


/ = -[ (o(d)— dx)=- f — (\ 6 ff(S)dS)dx = \ S ' ar(S)dS (3) 

Jr dx J r dx Jo Jo 

where 5, is the relative displacement at the crack tip. 

It has been shown that, for self-similar crack growth and negligible cohesive zones (or small-scale 
yielding), the /-integral can be expressed as the rate of decrease of potential energy with respect to the 
crack length and is equivalent to the energy release rate, G [31], [32]: 


J=-— = G 
da 


where II is the potential energy of the system. From (3) and (4): 


(4) 


G = [ S ‘ a(S)dS 

Jo 

Considering S,=S rnaA at crack extension the critical value of G, G c , is: 


(5) 


G c = J*”' <7(S) dS (6) 

When using Barenblatt’s [26] model in the simulation of delamination growth, the cohesive zone can 
still transfer load after delamination onset (point 3 in Figure 9), until the critical value of the energy 
release rate is attained (point 4 in Figure 9). 
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The concept of decohesion zones to simulate delamination growth in composites is usually 
implemented by means of decohesion elements connecting the individual plies of a composite laminate. 
These elements can model the discontinuity introduced by the growth of delaminations. Decohesion 
elements use a high penalty stiffness before delamination onset to prevent additional deformations (see 
point 1 in Figure 9). 

Decohesion elements can be divided into two main groups: continuous decohesion elements and point 
decohesion elements. Several types of continuous decohesion elements have been proposed, ranging from 
plane decohesion elements with zero thickness connecting solid elements [33-36]; plane decohesion 
elements with finite thickness connecting shell elements [37]; and line decohesion elements [38-40], Point 
decohesion elements are identical to spring elements connecting nodes [41, 42]. 

The concept of decohesion elements has been used in different types of problems: compression-after- 
impact [33, 35], stiffener-flange debonding [36], damage growth from discontinuous plies [40], 
diametrical compression of composite cylinders [37], The use of decohesion elements in the simulation of 
delamination growth in compression-after-impact problems is illustrated in Figures 10 and 11. Figure 10 
shows the deformed shape of a CFRP laminate containing a delamination. Figure 1 1 shows the initial 
delamination due to the impact load and its predicted propagation. Figure 1 1 indicates that decohesion 
elements can capture non-self-similar crack growth. 



Figure 10. 


Global buckling in a [0V90 2 4 ] CFRP laminate under compression-after-impact loading 
[35], 


9 








i /i/. # • * * 

*4 ■■#'■■#'*'■ 

^ f £ # f * # * 

1 4|* 4 * * * * * 

IfV****** 

|~#-*M>~4 t-#- 

* * * # 

||* 4 * * • 

.4W4-%4m*-- 



4*4 

4*4 

•4 •■*-•4 

* * * 

4*4 

4 ♦ # 

*hM* 

% * 4 
4 •* - ♦ 
* * * 
4 * 4 


-*-f 

* » 

■*■4 
* * 
•*•4 

* 4 

* * 
•*"4 

* 4 
♦••4 

* 4 
■-*-> 

* -4 

* 4 



Figure 1 1 . Delamination propagation for a (0 Q 4/90 e 4 ) CFRP laminated composite under compression- 
after-impact loading [35]. 


The following issues must be addressed in order to obtain accurate results for the simulation of 
delamination using interfacial decohesion elements. 


Constitutive Equations 

The need for an appropriate constitutive equation in the formulation of the decohesion element is 
fundamental for an accurate simulation of the interlaminar cracking process. In addition to Needleman's 
[28-30] interfacial behavior represented in Figure 8, other constitutive equations proposed are [43]: linear 
elastic-perfectly plastic, linear elastic-linear softening, linear elastic -progressive softening, linear elastic- 
regressive softening (Figure 12). 
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Perfectly plastic (pp) 

Linear softening (lin) 

~ Progressive softening (pro) 
Regressive softening (reg) 
Needleman (Ne) 



^lin ^Ne ^ reg 
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Figure 12. Constitutive strain softening equations [43], 


In order to simulate the cohesive zone process ahead of the crack tip represented in Figure 9, a linear 
elastic-linear softening behavior is usually implemented [33-39, 44-46], A high initial stiffness is used to 
hold the top and bottom faces of the decohesion element together in the linear elastic range. For pure 
Mode I, II or III loading, after the interfacial normal or shear tractions attain their respective interlaminar 
tensile or shear strengths, the stiffnesses are gradually reduced to zero. The area under the stress-relative 
displacement curves is the respective (Mode I, II or HI) fracture energy (see equation 6 and Figure 13). 
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( 7 ) 


j o ^<7 3 3(^3=Gc 
J^T u {S)dS 2 =G IIC (8) 

J 0 ^23 WS^Gnc (9) 

Therefore, the properties required to define the interfacial behavior are the initial stiffness (or 
penalty stiffness), K P , the corresponding fracture energies, G [C , Gu C and G n ic and the corresponding 
interlaminar tensile or shear strengths, <r 33 , f 13 and r 23 . 

Once a crack is unable to transfer any further load (point 5 in Figure 9), all the penalty stiffnesses 
revert to zero. However, it is necessary to avoid the interpenetration of the crack faces. Therefore, the 
contact problem is addressed by re-applying the normal penalty stiffness when interpenetration is 
detected. The contact condition is related to the constitutive relation proposed by Needleman [28] and 
shown in Figure 8. 

This interfacial constitutive equation, which is shown in Figure 13, can be implemented as follows: 

i) For Sj < Sbj , the constitutive equation is given by: 

~K P 0 0 " 

<r= 0 K p 0 d = DS (10) 

0 0 K P 

ii) For So,, < S < S max ,i , the constitutive equation is given by: 

g = (I- E)DS (11) 

where I is the identity matrix and E is a diagonal matrix defining the position of the integration point in 
the softening curve. 

iii) For S - Snmj , all the penalty stiffnesses revert to zero. Interpenetration is prevented by 
reapplying only the normal stiffness if crack closure is detected, meaning that frictional effects 
are neglected. 



Figure 13. Bi-linear constitutive equation. 


The process of reapplying the normal stiffness when interpenetration is detected is typical of solution 
procedures of contact problems using penalty methods in a constrained variational formulation. It is also 
clear from this formulation that the accuracy of the analysis depends on the penalty stiffness K P chosen 
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for the linear-elastic region of the constitutive equation. High values of K P will avoid interpenetration of 
the crack faces but can lead to numerical problems. 

Some authors have determined the value for the penalty stiffness as a function of the interface 
properties. Daudeville et al. [47], for instance, have considered the interface as a resin rich zone of small 
thickness, e„ and have proposed a penalty stiffnesses defined as: 


K 


i 

p 



2G 13 ^ w _ 2 G 23 

5 A p — 

e, e t 


( 12 ) 


where G 23 , Gu, and E 3 are the elastic moduli of the resin rich zone. Values that have been calculated for 
the penalty stiffness K P include, 5.7x 10 7 N/mm 3 [39], and 10 8 N/mm 3 [44], Other authors have shown that 
lower values such as 10 7 N/mm 3 [34] can adequately minimize the relative displacements at the interface 
while avoiding the potential for numerical errors related to computer precision. 

In order to fully define the interfacial behavior, the unloading response must be specified. Crisfield et 
al. [37, 38] and Daudeville [47] assumed that, with reversing strains, the material unloads directly toward 
the origin (see Figurel3). This procedure seems reasonable since the interfacial stiffness when reloading 
is lower than the original (undamaged) stiffness. Such a procedure simulates the effects of the previous 
damage mechanisms that occurred along the interface. Other authors [40] have proposed an unloading 
curve with a slope corresponding to Hooke’s law. Such a procedure, typically used in the formulation of 
plasticity problems, would lead to the use of the same penalty stiffness when reloading, and to permanent 
relative displacements along the interface when the load reverts to zero. 


Mixed Mode Delamination 

In structural applications of composites, delamination growth is more likely to occur under mixed- 
mode loading. Therefore, a general formulation for decohesion elements should deal with mixed-mode 
delamination growth problems. O’Brien et al. [12] have proposed using a cubic equation obtained from 
least square regression curve fit of the experimental data, as shown in Fig. 14. 
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Figure 14. Mixed mode delamination criterion for AS4/3501-6 [12] 
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The development of a criterion such as shown in Figure 14 requires a more extensive test program 
than is often available, hi the absence of such test data in the mixed-mode range, the analyst must rely on 
empirical expressions. A comprehensive study of failure criteria for mixed-mode delamination in brittle 
epoxy, tough epoxy and thermoplastic composites under the full mixed-mode range was performed by 
Reeder [48, 49]. Some of the criteria evaluated include: 


i) Linear Criterion: 


ii) Power Law Criterion: 


iii) Bilinear Criterion: 
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(13) 


(14) 


(15) 


Mixed-mode bending tests were used to measure the mixed-mode delamination toughness of the 
above composites, providing experimental data to assess the criteria. The linear model appears to be the 
most suited to predict failure of thermoplastic PEEK matrix composites because the results were very 
comparable to the more sophisticated criteria, while using two fewer independent variables. However, the 
linear criterion was found to be inaccurate for predicting the response of epoxy composites. The proposed 
bilinear criterion [48], taking into account the modification of failure mechanisms near the one-to-one 
ratio of GilGn, was found to yield the best results when simulating epoxy matrix composites. The power 
law criterion provides a conservative simulation of delamination in epoxy matrix composites. 

Chen et al. [38] have proposed the use of both the linear criterion and the power law criterion to 
predict total interfacial fracture under mixed-mode loading and obtained a reasonable agreement between 
predictions and experimental results in the simulation of fixed ratio mixed-mode (FRMM) tests when 
using a=P=4. Other authors [34] have used a criterion based on moving damage surfaces as a function of 
the normal and in-plane relative displacements. 

An issue that should also be addressed is the definition of a delamination onset criterion for multiaxial 
traction states. Under pure Mode I, II or III loading, the onset can be determined simply by comparing the 
traction components with their respective material allowables. However, under mixed-mode loading this 
approach must be recast, since delamination onset due to multi-axial traction states may occur before any 
of the traction components involved reach their respective allowables. 

Cui et al. [50] have shown that stress-based criteria can be used for uncracked specimens where there 
is no macroscopic singularity. Furthermore, the authors have shown that interaction between stress 
components is important in determining delamination initiation, since poor results were obtained by using 
the maximum stress criterion [50], This observation is confirmed by Mohammadi et al. [51], who state 
that it is widely accepted that the quadratic failure criterion based on interlaminar stresses can be properly 
used to predict the initiation of delamination. 


Element Integration 

Special care must be taken in the integration scheme adopted for decohesion finite elements as 
oscillations of the traction field may occur [34, 39, 44, 45, 52], Schellekens and de Borst [44, 45] have 
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demonstrated using eigenmode analysis of the element stiffness matrices that Gaussian integration can 
cause undesired spurious oscillations of the traction field when large stress gradients are present over an 
decohesion element. With increasing mesh refinement, the element performance improves due to the 
decreasing stress gradients. 

The authors concluded that for line decohesion elements and linear plane decohesion elements the 
performance can be improved by using either a nodal lumping scheme, Newton-Cotes, or Lobatto 
integration scheme. It was shown that the eigenmodes of elements integrated with such techniques do not 
show coupling between displacements at adjacent nodes [45]. 

For quadratic decohesion elements, Newton-Cotes and Lobatto integration schemes produce smooth 
traction profiles when the displacement field over the element varies in only one direction. Oscillatory 
traction profiles may occur in other cases [52], However, this adverse effect is only significant for the 
central integration point, and it is less pronounced than the oscillations occurring when Gauss integration 
is used. If necessary, more refined meshes or linear elements can be used to reduce the magnitude of the 
stress oscillation. 

Another relevant issue related with the integration of decohesion elements is the use of full integration 
schemes. Analyses of problems involving crack propagation and strain-softening behavior have shown 
that the use of full integration was superior to the use of reduced integration schemes [53]. However, 
Alfano and Crisfield [54] have shown that for fully integrated linear 4-node decohesion elements, 
increasing the number of Simpson integration points from 2 to 20 results in an increase of spurious 
oscillations in the load-displacement curve and, consequently, to less a robust solution algorithm. 


Non-linear Solution Procedures 

The softening nature of the decohesion element constitutive equation causes convergence difficulties 
in the solution of the Finite Element Analysis. Crisfield et al. [43] found that when using the Newton- 
Raphson method, under load (with the arc-length method) or displacement control, the iterative solutions 
often oscillated when a positive slope of the total potential energy was found, and therefore failed to 
converge. In order to obtain convergence, a line search’ procedure with a negative step length was 
proposed. 

In a recent work, Mi et al. [39] proposed the use of a modified cylindrical arc-length method [55] to 
obtain converged solutions. This technique chooses between the alternative roots to the arc-length 
constraint, selecting the one involving the minimum residual [55], However, even with this new 
technique, the solution sometimes enters a cycle and oscillates between two points. The oscillation must 
be detected so that the increment size can be reduced. It was also found that when coarse meshes are used, 
a ’snap-back’ type of behavior might occur [39], Therefore, in order to obtain a relatively smooth solution, 
the mesh should be fine enough to include at least two decohesion elements in the cohesive zone located 
at the crack tip. This means that for Mode I loading and for a given G/c, the lower the interlaminar tensile 
strength <7 33 , the larger will be the cohesive zone, and hence the smoother the solution. This 'root 
selection’ method has also been successfully used by Wisheart and Richardson [46] in the simulation of 
DCB tests in composite laminates. 

In the previously mentioned work, Chen et al. [38] have used two different techniques to obtain 
converged solutions: displacement control in conjunction with the line searches in the commercial finite 
element code ABAQUS and the modified cylindrical arc-length method [55] in the LUSAS finite element 
code. It was found that it is significantly easier to obtain convergence when using linear decohesion 
elements. A LUSAS model using quadratic elements failed to converge. Similarly, an ABAQLTS model 
with the same interlaminar tensile strength and linear elements also failed to converge. The convergence 
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difficulties were overcome by reducing the interlaminar tensile strength to 40% of the original value 
whilst maintaining a constant fracture energy. The improvement in convergence is due to the previously 
mentioned effect of increasing the cohesive zone ahead of the crack tip. This modification had only an 
effect on the behavior around the ultimate load of the structure. As expected, when decreasing the 
interla mi nar tensile strength, the predicted peak load also decreased. 

The arc-length method proposed by Riks [56] to pass limit points in non-linear problems uses a 
constraint equation to relate the incremental load factor to the norm of the incremental displacement 
vector. However, Schellekens [51] recognized that the arc-length method has difficulty converging in 
problems such as the simulation of delamination using decohesion elements where failure is highly 
localized. Therefore, the authors suggested that the displacement norm should be determined considering 
only the dominant degree of freedom, as the definition of a constraint equation based on the global 
displacement vector is less effective. A nodal crack opening displacement (COD) control was proposed 
[51] where the load factor is calculated from a reduced displacement vector that is a function of the COD. 

Reddy Jr. [37] used explicit time integration to overcome the difficulties of having to deal with 
material softening and contact conditions when using decohesion elements. Although explicit 
formulations are typically used to solve transient dynamic problems, they were used to address a problem 
involving quasi-static loading. The explicit time integration methods are conditionally stable because the 
minimum time step used for the explicit time integration of the governing equations depends on the 
highest eigenvalue in the mesh. Estimates for the frequency of the vibrational modes that influence the 
choice of the stable time increment were performed and it was found that the frequency depends on the 
penalty parameter K P . However, this analysis is valid before delamination onset, as it is expected that 
when the element is softening the minimum time step used for the explicit time integration is constantly 
being modified. 

Ladeveze et al. [57-59] have proposed a 'large time increment method' to solve the nonlinear problem. 
The procedure is based on a single global iterative procedure on the whole loading history between linear 
steps and non-linear steps. 

Discussion 

In order to improve the generality and the accuracy of the formulations described, a number of 
developments are required. These may include: 

- A criterion able to predict delamination onset under multi-axial stresses. The use of mixed-mode 
criteria to predict total element decohesion is proposed. 

- Decohesion element constitutive equation: Experimental results have shown different damage 
mechanisms under pure Mode I and Mode II loading. However, most of the approaches assume a 
bi-linear constitutive equation for both loading modes. An investigation of the appropriate 
constitutive equations for different loading modes, including unloading behavior, is required. 

- The development of decohesion elements compatible with plate/shell elements would allow the 
effective analysis of large-scale structures (e.g. debonding in skin-stiffener). 

- Usually, the fracture process of a composite structure occurs by a combination of intralaminar (e.g. 
matrix transverse cracking, fiber fracture) and interlaminar damage (Figure 1). Therefore, a model 
incorporating the simulation of the different damage mechanisms would provide a general tool for 
the failure analysis of a composite structure. 

- Point-to-surface contact algorithms such as the one proposed in [34] should be developed to avoid 
interpenetration of the crack faces after complete decohesion, especially when large tangential 
relative displacements occur. 
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Concluding Remarks 


The VCCT is a valuable technique for the determination of both delamination growth and the 
contributions of different loading modes to the global energy release rate. Important information can be 
obtained without having computationally expensive models, using only interlaminar fracture toughness 
for del ami nation propagation prediction, and no stress-based allowable. However, the VCCT must use a 
pre-defined crack and for some structures/loads the precise location of crack initiation might not be 
obvious. Furthermore, complex moving-mesh techniques are required to simulate delamination growth. 

The use of decohesion elements can overcome some of the above difficulties. Using decohesion 
elements, both onset and propagation of delamination can be simulated without previous knowledge of 
crack location and propagation direction. 

The material properties required to define the element constitutive equation are the interlaminar 
fracture toughnesses and the corresponding strengths. The fact that the strengths are required is a 
disadvantage when compared with the VCCT, which only requires interlaminar fracture toughnesses. 

It is also clear that Newton-Cotes integration schemes should be used in the definition of the stiffness 
matrix of decohesion elements because Gaussian integration may lead to oscillations of the stress field. 
Due to the localized softening associated with delamination growth, the solution procedures for the 
nonlinear problem requires special techniques such as the COD control. Such techniques can enhance the 
numerical stability and, consequently, improve the convergence rate. 
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